Structure and dielectric properties of polar fluids with extended dipoles: 

results from numerical simulations 



cn 

O 

o 

> 
O 

:^ 

oo 
(N 



V. Balleneggei0 and J. -P. Hansen 

Department of Chemistry 
University of Cambridge 
Cambridge CB2 lEW (UK) 
(Dated: February 2, 2008) 

The strengths and short-comings of the point-dipole model for polar fluids of spherical molecules 
are illustrated by considering the physically more relevant case of extended dipoles formed by two 
opposite charges ±g separated by a distance d (dipole moment ^ = qd). Extensive Molecular 
Dynamics simulations on a high density dipolar fluid are used to analyse the dependence of the 
pair structure, dielectric constant e and dynamics as a function of the ratio d/a {a is the molecular 
diameter), for a fixed dipole moment /x. The point dipole model is found to agree well with the 
extended dipole model up to d/a ~ 0.3. Beyond that ratio, e shows a non-trivial variation with d/a. 
When d/a > 0.6, a transition is observed towards a hexagonal columnar phase; the corresponding 
value of the dipole moment, fi^ /a^kT = 3, is found to be substantially lower than the value of the 
point dipole required to drive a similar transition. 
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I. INTRODUCTION 

Highly polar fluids are particularly important in many 
areas of Physical Chemistry, Chemical Engineering and 
Biology, because of their role as solvents leading to elec- 
trolyte and polyelectrolyte dissociation. Water is of 
course the most important among polar liquids, but be- 
cause of its complex behaviour, primarily linked to the 
formation of hydrogen-bond networks, much theoretical 
work has focussed on simpler models involving spheri- 
cal molecules with point dipoles. The best-known and 
widely studied examples are dipolar hard spheres, and 
the Stockmayer model (dipolar + Lennard-Jones inter- 
actions). A long-standing problem, going back to the 
classic work of Onsager and Kirkwood 0, is to re- 
late the dielectric response of a polar fluid to molecular 
dipole fluctuations and correlations Q. Subtle concep- 
tual and numerical problems arise in Molecular Dynamics 
or Monte Carlo simulations of finite samples of polar flu- 
ids, which are linked to the infinite range of the dipolar 
interactions, so that boundary conditions must be treated 
adequately. These issues were resolved in the early eight- 
ies, both for the reaction field and the Ewald summation 
implementations of boundary-conditions 0, 0, 0- I-*^" 
spite this theoretical progress, accurate estimates of the 
dielectric permittivity of simple polar fluids by numerical 
simulation remains a very challenging task, because large 
fluctuations of the total dipole moment of the sample oc- 
cur on a relatively long timescale (of the order of 10 ps), 
leading to a very slow convergence rate for the dielectric 
constant [1,0] (see also Set. HI). 

More recently, it was realized that simple dipolar liq- 
uids can form a ferroelectric nematic phase for sufficiently 
large dipole moments [Tol HTl IT^ . This transition is in- 
timately related to the formation of chains of dipoles 



aligned head-to-tail, which prevent the formation of a 
proper liquid phase in the Stockmayer model if the dis- 
persive energy is below a certain threshold . 

However point dipoles represent a limiting situation, 
never achieved in real polar molecules, which are char- 
acterized by extended charge distributions linked to 
electronic charge transfer from electron donors to elec- 
tron acceptor atoms. In simple heteronuclear diatomic 
molecules like CO or HF, this situation can be modelled 
by assigning fractional charges of opposite sign to sites 
which are separated by a distance d, typically of the order 
of 0.1 nm [iSj. Such, or more complicated situations in- 
volving more than two atoms, can be mimicked by adding 
higher order point multipoles to a point dipole Jl5j|, but 
such an expansion will require more and more high order 
multipoles as two molecules approach each other. 

In this paper, we present a systematic investigation 
of the structure, dielectric response and phase behaviour 
of a simple model involving spherical molecules carry- 
ing extended (rather than point) dipoles resulting from 
opposite charges ±g, each displaced symmetrically by a 
distance d/2 from the centre of the molecule. We study 
how the properties of the polar liquid change when d is 
increased from zero, varying q simultaneously so that the 
dipole moment = qd remains constant. Although po- 
lar molecules are never spherical, the model investigated 
in this paper, which focusses on the electrostatic rather 
than steric interactions, is the simplest "natural" exten- 
sion of the dipolar sphere model towards a more realistic 
representation of highly polar fluids. Some studies on the 
structure of similar models with extended dipoles have 
been published previously, but without an investigation 
of their bulk dielectric properties [l^ IT^ . 



II. THE MODEL AND SIMULATION DETAILS 
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We consider a polar fluid made up of soft spheres with 
two embedded point charges ztq located at ±d/2 from 
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the centre of the molecule (see Fig. The distance |d| 
is assumed fixed, so the molecule is not polarizable and 
carries a permanent dipole moment fi = qd. 




FIG. 1: A polar molecule with an extended dipole moment 



Placing the origin at the centre of the sphere, the 
multipole moments qim = / ^jml^*, '/')^'/o(i")d'^i"j where 
p{r) = q6(d/2) — q6{—d/2) is the molecular charge dis- 
tribution nal, are 



qir. 




if I odd and m 
otherwise. 



(1) 



The next non-vanishing moment after the dipole is thus 
the octopole, since the quadrupole moment vanishes by 
symmetry for this choice of origin. 

The interaction energy between two molecules is given 
by the sum of a soft sphere repulsion Wss(r) — 4utT^^/|r|^^ 
and the four Coulombic energies due to the point charges. 
On Fig. 13 the electrostatic energy at contact is com- 
pared to a truncated multipolar expansion containing the 
dipole-dipole and dipole-octopole interactions. The con- 
figuration of lowest energy occurs when the molecular 
dipoles are aligned head-to-tail {9 = 02 = 0). This mini- 
mum energy is lower for extended than for point dipoles. 



exact Coulombic energy 
dipolar interaction 
dipolar + dipole-octopole int. 




FIG. 2: Electrostatic interaction energy of two molecules at 
contact (|r| = a) for d = a/2. 



The polar fluids we have considered can be character- 
ized by the following dimensionless (36i] parameters: 



• reduced density: 

• reduced temperature: 

• reduced dipole moment: 



pa^ = 0.8 
= kT/u = 1.35 



/i* = y/p'^/a^u ^ 2.0 
• reduced moment of inertia: /* = I/ma^ = 0.II7. 
This corresponds to the same thermodynamic state point 
that has been extensively studied by Kusalik in the case 
of point dipoles See '2^ for an investigation of the 

phase diagram of the dipolar soft-sphere system. Equi- 
librium quantities, such as the dielectric constant and 
distribution functions, are independent of the moment of 
inertia. 

In all calculations, we employed periodic boundary 
conditions. We chose the spherical geometry, i.e. the 
periodic replications of the basic cubic simulation cell 
form an infinite sphere, which is itself embedded in a in- 
finite region of dielectric constant e'. In this case, the 
Hamiltonian of the system is 



N 



i<3 = l 



27rM^ 



(2e' + 1)L3 



(2) 



where L is the side of the box, M — J2i Qi^i is the total 
dipole moment, and 



crfc(K|r + nL\ 



-y- 

n=iO 



nL 



exp 



27ri 

] 

L 



(3) 



The last term in jSJ accounts for the work done against 
the depolarizing field due to the surface charges induced 
on the spherical boundary. This term vanishes only for 
metallic boundary conditions (e' = oo). The Ewald sums 
in 'l'(r) were evaluated using the Smooth Particle Mesh 
Ewald method p3| (Ewald coefficient k = 3.4705 nm~\ 
grid size 32x32x32, interpolation order 6). The interac- 
tions were truncated beyond 0.9 nm, both for the real 
space Ewald sum and for the soft-sphere repulsions. 

Molecular Dynamics simulations were carried out us- 
ing the simulation package gromacs "22]. The equations 
of motion were integrated using the so-called leap-frog al- 
gorithm with a reduced time step of dt* ~ dtj \J ma'^ ju = 
0.00235. The temperature was kept constant using a 
Berendsen thermostat. Equilibration periods lasted at 
least 100 ps (~118 reduced time units), and were followed 
by data producing runs of 8 ns or more. The number of 
molecules was 512 in calculations of the dielectric con- 
stant (Set. Ill), and 5555 in calculations of correlation 
functions (Set. IV). 



III. DIPOLE FLUCTUATIONS AND 
DIELECTRIC CONSTANT 

The dielectric constant of a homogeneous and isotropic 
fiuid can be calculated from the fluctuation formula 



(e- l)(2e^ + l) _ 4^ (M^) 



2e' 



31^ kT 



(4) 
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diffusion constant was calculated from Einstein's relation 
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FIG. 3: Dielectric constant of a dipolar soft sphere fluid as a 
function of dipole extension (p* — 0.8, T* = 1.35, /x* = 2). 



which holds for a macroscopic spherical sample of vol- 
ume V surrounded by a medium of dielectric constant e'. 
The results obtained for the dielectric constant are inde- 
pendent of the choice of e', provided the boundary term 
in eq. |5J) is properly taken into account. We employed 
metallic boundary conditions, because they are known to 
produce less uncertainties in estimates of e than bound- 
ary conditions with finite e' |HI13 (see also below). The 
fluctuation formula reduces in this case to 



e=l + 3y{g), g = 



TV' 



(5) 



where the dimensionless parameter y = 47r/3p/i^/9 equals 
3.309 at the state point under consideration. 



d/a 



TM [ps] [ps] D [10~^ cm7s] 



0.02 99.6(±1.4) 


2.21 


0.50 


11.7 


0.50 


0.3 98.4(±1.5) 


2.64 


0.54 


11.6 


0.53 


0.4 94.0(±1.5) 


2.88 


0.63 


11.5 


0.53 


0.5 92.4(±1.7) 


4.14 


0.88 


10.6 


0.49 


0.6 102.3(±3.2) 


11.44 


2.01 


8.7 


0.31 


0.61 104.7(±3.6) 


13.97 


2.36 


8.5 


0.27 


0.62 100.8(±3.5) 


14.81 


2.79 


7.9 


0.22 



TABLE I: Influence of dipole elongation on some properties 
of a polar soft sphere fluid. 

We show in Fig. O the variation of the dielectric con- 
stant with the extension d of the dipole, at fixed dipole 
moment /z* = 2. Table HI lists the actual values of e, 
together with some other quantities of interest, namely 
the diffusion constant D, the dielectric relaxation times 
tm and r^, and the reduced pressure p* = pa^/u. The 



(|r,(t)-r,(0)|2) =6I?i, t 



(6) 



The relaxation times tm and were determined from 
the autocorrelation functions C]v[(i) = • M(0)) / 

(AP) and C^it) (see Fig. EJ. For t > 0.3ps, CmW 
exhibits an exponential decay exp(— t/rM) typical of a 
Debye dielectric. The relaxation oiCfj,{t) is not well fitted 
by a single exponential, and the corresponding relaxation 
time was estimated from the integral of C^i (t) . 




FIG. 4: Autocorrelation functions Cmit) and C^{t), for the 
value d — a/2. The inset shows a logarithmic plot, confirming 
the exponential behaviour of Cm{t). 

For almost point dipoles {d* = d/a = 0.02), our value 
for the dielectric constant is in good agreement with the 
result epoint = 98 ± 2 reported by Kusalik et al. for 
the ideal dipolar soft sphere fluid (25|]- Our data shows 
that when d* increases, the dielectric constant decreases 
and reaches a minimum about 6% lower than epoint at 
d* ~ 0.55. When d* is further increased, the dielectric 
constant increases rapidly above Cpoint, up to the criti- 
cal distance d* ~ 0.63. At this critical distance, a phase 
transition occurs from an isotropic fluid to an orienta- 
tionaly ordered "liquid crystal" phase (see SctlV|l. 

The simulations show that the point dipole model gives 
a reliable estimate of the dielectric constant over a very 
wide range of extensions d, namely up to the point where 
the system undergoes a phase transition. The weak sen- 
sitivity of the dielectric constant on the extension of 
the dipole, which contrasts with the large sensitivity ob- 
served in water models 0], may be due to the absence 
of a quadrupole moment in our molecules. 

It is clear from Table that the dynamics of the fluid 
slows down when d is increased: the diffusion coeffi- 
cient D drops and the relaxation times increase. This 
slowdown is due to the formation of head-to-tail dipo- 
lar chains in the system. Their entanglement make 
these chains less mobile than individual molecules in the 
present high density regime. 
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FIG. 5: Convergence of e with simulation time, for dipole 
elongations d* = d/a = 0, 0.3, 0.4, 0.5 and 0.6. 



Long runs were needed to obtain even moderate ac- 
curacy (about 2%) in the estimated dielectric constants. 
Fig. 13 shows the running estimate of e as a function of 
simulation time. The slow convergence, especially for 
large elongations of the dipole, can be traced back to the 
long relaxation times tm , as shown by the following error 
analysis. 

By definition, the probability distribution of the sam- 
ple having a total dipole moment of magnitude M and 
arbitrary orientation is given by 



P{M) oc AnM'^e 



-l3F{M) 



(7) 



where F{M) is the free energy of the system. From 
macroscopic electrostatics, the energy of a spherical di- 
electric sample, of dielectric constant e and carrying a 
uniform polarization M/V, is 



U{M) 



2e' 



V (e-l)(2e' + l)' 



(8) 



where e' is the dielectric constant of the surrounding 
mcchum. Using F{M) ~ F(0) + U{M), we find from 
and that the fluctuations of 5 = AP/{N^i^) are 
distributed according to 



9y 



2e' 



2 (e- l)(2e' + l)^ 



(9) 



where the normalization constant is yl = 2^ k-^/tt. The 
mean of this distribution is (g) — 3/ (2k), in agreement 
with the fluctuation formula Though the distribu- 
tion ((nj is approximate, since we neglected changes in 
entropy and eq. |(SJ) is valid only in the linear regime, it 
gives a good description of fluctuations of the total dipole 
moment observed in simulations of polar fluids |2^ . 

The dielectric relaxation time tm gives a time scale 
for two measurements of M'^ to be independent. In a 



simulation of total duration t, the distribution ij^l is 
thus sampled n ~ I/tm times. After n such indepen- 
dent measurements, the standard deviation in the aver- 
age °f 9 factor is cTg^n = Cg/v^ where 
— ((5 ~ (5))^) = 3/(2k^) is the variance of the dis- 
tribution The expected relative uncertainty in the g 
factor is therefore 



Ha) 



(5) 



Solving for e, one has 

3y (g) {2e' + 1) 
2e' + 1 - 3y (5) ' 



1 = 



(10) 



(11) 



By the rules of propagation of errors, the relative uncer- 
tainty in the dielectric constant minus one is thus 



2e' 



2e' + l 



3~' 



(12) 



The error bars in Fig. O were determined from this for- 
mula, and are in agreement with the fluctuations ob- 
served in Fig.|31 In a Debye dielectric, the relaxation time 
Tm is related to the Debye relaxation time td (which is 
independent of boundary conditions) by |2^ 



2e' + 1 



TM 



2e' 



-TD- 



(13) 



Inserting H13|) into (|12|l , we see that larger values of e' will 
lead to smaller uncertainties in the dielectric constant. 
This explains the faster convergence of e observed when 
using metallic boundary conditions p, '23] . 

According to the present analysis, the slow convergence 
of e, as determined from the fluctuation formula, is due to 
the large value of the Debye dielectric relaxation time 
and the rather broad distribution P{g). Moreover, the 
uncertainties in e are independent of system size, as long 
as it is macroscopic. In large systems, it may therefore be 
more favourable to determine e from correlation functions 
rather than from the fluctuation formula. 



IV. STRUCTURE 
A. The pair distribution function 

The pair distribution function h{\, 2) = /i(r, /Xj^, 1^2) of 
the infinite system can be expanded in rotational invari- 
ants (28|: 

/i(l, 2) = /i"°"(r)+/(,"°(r)$""(l, 2)+/i"2(^)^ii2(^^ 3) + . . 
where 

$"«(l,2)=./ii-/i2 (14) 
$"2(1, 2) = 3(Ai • f)(/i2 • f) - Ai • A2- (15) 
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The (I)'i'2i'g form on orthogonal basis for the angular de- 
pendence of h{l, 2). The first projections are 



'(r) = 3(Ml,2)$"°(l,2))^^_ 
^(r) = ^(Ml,2)$"^(l,2))^^^ 



M2 
M2 



(16) 
(17) 

(18) 



where {■ ■ ■) ^ = J ■ ■ ■ dil^/47r denotes an unweighted an- 
gular average over the orientations of /x. 

Plots of /i*'°''(r) and h^^'^{r) are shown on Fig. for 
three elongations d of the dipole. As d is increased, the 
stronger multipolar moments carried by the molecules 
lead to a slight reduction of the fluid structure as 
measured by the centre-to-centre distribution (?(r), but 
more orientational order, as measured by the projections 
h^^^{r) and h^^'~'{r) (the latter projection, not shown in 
the figure, closely resembles h-^^'^{r)). 




size of the volume element is properly taken into account 
in the normalisation |^ ^| . The reason for this sharp 
drop is simple: the Ewald potential deviates strongly 
from the Coulomb potential at distances greater than L/2 
(it vanishes in particular at a finite distance between L/2 
and \JiLj2 depending on the position of the charges). 
The rapid decrease in r^h}^'^(r) reflects thus an artiflcial 
flnite size effect arising from the Ewald summations of 
the periodically replicated system ||33] • In a system with 
long-range forces, great care must therefore be exercised 
in the interpretation of correlation functions at distance 
larger than half the box length. 
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FIG. 7: Convergence of r^h}^'^(r) at large distances towards 
the limit l)19|l . Data from a 6ns long simulation of a system 
of 5555 molecules (p* = 0.8, T* = 1.35, = 2, d* = 0.5). 



FIG. 6: Projections h°'^'-\r) and ft"^(r) of the pair correlation 
function for three values of d* = dja. 

The projection h^^'^ir) is related to the dielectric con- 
stant of the fluid by the formula 

lim r3/iii2(^) ^ 1 (19) 

first derived by Nicnhuis and Dcutch 27]. A 512- 
molecule system with a half box size of L/2a = 4.3 is 
too small to reach the asymptotic limit (|19|l . The results 
for the correlation function shown below in Figs. IT HlOl 
were hence obtained using a larger system (simulation of 
5555 molecules during 6 ns) under the same conditions 
(p* = 0.8, T* = 1.35, fi* = 2, d = a/2, e' = oo). Now 
L/2a — 9.55, and Fig.[7|shows that r^h^^^{r) does reach 
the asymptotic value 1191) at a distance r ~ 7a, as in the 
case of point dipoles [l3- Contrary to the fluctuation 
formula, eq. H19|l gives estimates of e that become more 
accurate when the size of the system is increased. 

Kusalik made the observation that r^/i^^^(r) drops 
sharply for r greater than L/2, even when the reduced 
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FIG. 8: The function r-^(/i"''(r) - /i^^°(oo)) and its integral 
(dashed line) . Same system as in Fig. Q 

The projection h^^°{r) is also related to the dielectric 
constant, since the fluctuation formula I^J can be written 



6 



in the Kirkwood form l38 



{e~l){2e' + l) 
2e' + e 



3^1 



IT 



/i,V"(r)rMr . (20) 



Since the l.h.s. of eq. H2(J|I depends on e', the projection 
h^^^{r) must also be sensitive to this boundary condition, 
whence the introduction of a subscript e'. Fig. |S1 shows a 
plot of r'^h]^{r), and the integral of this function, for the 
same system as in Fig. The correlations extend up to 
r ~ 7(T, just as in the case of h^^'^{r). This distance cor- 
responds to the scale beyond which the fluid behaves as 
a continuum dielectric and obeys the equations of macro- 
scopic electrostatics. 

Though formula H2U|) is equivalent to l@J, it is worth- 
while to discuss how the pair correlation function de- 
pends on the dielectric constant of the external medium. 
This problem has been addressed in Ref. [23| (see also the 
perturbation theory presented in Ref. |^); here, we hope 
to give a clear and concise answer to the above question, 
using simple physical arguments to justify the formulae. 

In a spherical sample, h\}^{r) is in fact the only pro- 
jection, among all /I'^'^^'s, to be strongly affected by the 
boundary condition e'. This is due to the surface term in 
the Hamiltonian Q , which corresponds to an interaction 
energy between two molecules 



47r Hi ■ fi2 

2e' + l V 



(21) 



which has the angular dependence of (f>^^"(l, 2). 

As the interaction 121|l is independent of the distance 
between the molecules, h\}^{r) does not decay in general 
to zero at infinity, but rather to an e'-dependent constant. 
We will prove below that this constant is 



lim Kr{r) 



2 (6-1) 



V 3pye 2e' + e 



(22) 



in the spherical geometry. The constant l|22|) vanishes 
only when using the boundary condition e' — e, which 
mimicks an infinite sample, or in the thermodynamic 
limit V oo (which is never reached in simulations). 

The fact that /i^,^°(r) contains the 0{1/V) constant 
contribution (|22|l at large distances ensures that the Kirk- 
wood formula (|2()|l gives consistent results for different 
boundary conditions. Indeed, when hl}'^{r) is integrated 
over the volume V = AttR^ /3 of a large sample, as in the 
r.h.s of eq. (|20l) . the constant H22() gives a finite contribu- 
tion to the integral that is precisely what is required to 
match the e'-dependence of the l.h.s. of the equation. In 
other words, eqs. (|21|l and (|22|) imply that 



47r 



,110 



(r)r^dr-fy/i^V"(oo), 



when the samples are large enough (i.e. in the limit R 
oo). When this identity is inserted into the Kirkwood 
formula (|2()|l . it is immediately clear that the predicted 
dielectric constant is independent of e'. 



In order to prove eq. H22() with simple physical argu- 
ments, we need to recall two basic results from the sta- 
tistical mechanics of polar liquids. The first result is 
the expression of the density p{r, fi) of molecules at r 
with orientation /x in a polarised sample permeated by a 
macroscopic field E(r): 



1 



P(r,/.) = ^(1 



32/ 



Pfi-E{r)] +OiE') (23) 



{y is defined after eq. ©). This formula is consistent 
with the constitutive relation P(r) = (e — l)/(47r)E(r) 
of macroscopic electrostatics, since the average polar- 
isation density in the fluid is by definition P(r) = 

The second result we need is the expression of the ef- 
fective dipole moment /i,®* of a polar molecule held fixed 
in a polar liquid (/.t°^ is defined as fi, the dipole moment 
of the fixed molecule, plus the total dipole moment of the 
screening cloud around /x). One may first think naively 
that fi'^^ = fi/e: the fluid would screen the dipole accord- 
ing to its dielectric constant. But this would be treating 
the polar fluid as a dielectric continuum everywhere, in- 
cluding in the interior of the flxed molecule, which is 
obviously wrong. An exact statistical mechanical calcu- 
lation shows that the right answer is j30| 



e- 1 
3ye 



(24) 



(The expression (e — l)/3ye can itself be interpreted as 
being composed of a factor (e — l)/32/ arising from lo- 
cal correlations around fi, times the expected factor 1/e 
due to screening by distant molecules). With these two 
results in mind, we can now understand easily formulae 
(EH, (EB and mi). 

The result (|19() for the large distance behaviour of the 
pair correlation function h(l, 2) can be seen as a straight- 
forward consequence of the screening effect H24f) . By 
definition of the distribution functions, the density of 
molecules at r2 with orientation fj,2, when a molecule 
is known to be located at 1 = (ri,/i]^) is 



P(2|l) = 



P(2,l) 
P(l) 



Mi,2)). 



(25) 



From l|24|l . the electric field due to the fixed molecule fii 
is equivalent, at large distances ri2 = r2 — ri, to that 
of a renormalized dipole moment . This dipolar field 
— V2(Ati^ • Vi)l/|ri2| is locally uniform and weak, so we 
can apply the linear response result H23|l . Using H23(l and 
we find that 



P(2|l), 7 ^(l^^^-4^Pv^,^{l,2) 



(26) 



where Wdip( 1,2) = (/x^ • Vi)(/22 • V2)l/|ri2| is the dipolar 
potential, and we assume an infinitely extended sample. 
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Comparing (|25|l and (|26|l . we obtain the asymptotic be- 
haviour of the pair correlation function: 



h{l,2) 



(-/3t)dip(l,2)), |ri2|->(». (27) 



The result follows then upon inserting ((23 into ((TH|l . 

Formula |(23) can be interpreted in a similar way, 
namely as arising from the interaction of a molecule with 
the reaction field produced by the screened dipole mo- 
ment of another molecule. We recall from macroscopic 
electrostatics that a dipole /ij^ at the centre of a spherical 
sample of radius R and dielectric constant e, surrounded 
by a dielectric medium e', produces a uniform reaction 
field 



Er'(/^i) 



e 2e' 



i?3 



(28) 



inside the sample, because of the surface charge density 
induced at the dielectric discontinuity In a finite 

spherical sample, a molecule at a position Y2 far enough 
from - so that it does not disturb the screening cloud 
around it - will interact therefore not only with the dipo- 
lar field of /x^*^, as in (|26|l . but also with the reaction field 
E_r(/^i*^) produced by the screened dipole moment of this 
molecule. From (|23|l , a term 



47r 



1 



3y 



(29) 



must hence be added to (|26|) in a finite spherical sample. 
The pair correlation at large distances, eq. H27() . includes 
then the additional contribution 



4^ 2(£- 1)^ e'-e 
W 9y2e 2e' + e 



/3/X2 • ^ll 



(30) 



where we used and V = 47ri?^/3. Formula if^ 

follows now from projecting (|30() onto $^^"(1, 2) =(1^ (12 
according to l(T7|l . 

We conclude this discussion by noting that the interac- 
tion energy H21|) used in the simulations, or equivalently 
the surface term in the Hamiltonian |(2Jl, can also be in- 
terpreted in terms of a reaction field effect. Indeed, each 
polar molecule in the sample will create a reaction field, 
acting on itself and on all other molecules, that is given 
by eq. (^5)1 with e = 1 . Since the Ewald sums ^ give the 
electrostatic energy between the molecules in the case of 
a sample surrounded by a metal (e' — cxd), the correc- 
tion to this energy to be used in the Hamiltonian of a 
spherical system with boundary condition e' is 



1 



N 



^(-M.)- (e1^'^1(m,)-E[^^°°1(m,) 



in agreement with |(2Jl. 



27rM^ 
e' + l 



, (31) 



B. Site-site correlations 

Contrary to the point dipolar fluid model, the present 
model with extended dipoles has well defined site-site 

distribution functions /i++(r) = /i__(r) and (r) [s^ . 

From these, we get the charge-charge correlation function 

S{r) = S'intra(?') + S'i„tcr(?'), whcrC 



(32) 



describes the inter-molecular correlations, while 

.(r) = 2ep5iv) ~ 2gV^^fe^ (33) 



^ ' 47rd2 

is the intra-molccular correlation function for a molecule 
with a rigid dipole of extension d. The charge-charge 
correlation is of sp ecial interest, because it satisfies the 
two sum rules [33 : 



Neutrality: j S{r) d^r = (34) 
Stillinger-Lovett: 1 = 1 + ^ j d?vY^ S{r). (35) 




FIG. 9: Site-site distribution functions. Same system as in 
Fig. 13 



The site-site correlation functions and S{r) are shown 
on Fig. for d — a/2. The charge neutrality sum rule is 
found to be accurately satisfied: 



{h++{r) - h+^{r))r^<:\r ~ 8.8 • 10"^ (36) 



The Stillinger-Lovett sum rule allows in principle the de- 
termination of the dielectric constant from S'(r), but this 
route is not practicable in a computer simulation, be- 
cause of the unfavourable ratio (1 — e)/e which saturates 
for large e, and also because it is difficult to determine 
the second moment of 5'intor('') accurately. Fig. 1 1 01 shows 
however that eq. (|35|l is satisfied within the uncertainties 
of our data. 
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0.02 - 



0.01 



-0.01 



-0.02 - 



I I , I r- 

4^/;S'intcr(r')r-'Mr-' 




d/a 




Pi 


P2 




e± 


0.62 


0.80 


0.08 


0.07 


e = 


103.8(±4) 


0.63 


0.80 


0.09 


0.08 


e = 


112.6(±5) 


0.64 


0.95 


0.97 


0.91 


1.16 


1.46 


0.65 


0.98 


0.66 


0.91 


1.28 


1.56 


0.66 


1.03 


0.98 


0.94 


1.02 


1.43 



TABLE II: Constant pressure simulations of the dipolar soft- 
sphere system at p* = 0.22, T* — 1.35 and fi* — 2. Data from 
8 ns long simulations of 512 molecules, collected after an equi- 
libration period that lasted up to 10 ns. For d > 0.64(t, the 
system is a ferro-electric liquid crystal with columnar order. 
Uncertainties in ey and e± are about ±0.01 and ±0.02 respec- 
tively. 



FIG. 10: Integration of the second moment of Sinter (»")• 



V. ORIENTATIONAL ORDER 

When the molecular dipole has an extension greater 
than d ~ 0.64a, the simulations show spontaneous for- 
mation of orientationaly ordered phases, starting from 
random initial configurations. At the state point under 
consideration {p* = 0.8, T* = 1.35), we observed phase 
coexistence between a dense liquid crystal and a very di- 
lute gas. In order to deal with pure phases, we performed 
simulations at constant pressure {p* = pa^/u = 0.22), 
rather than constant volume, for d > 0.62(7. 

The occurrence of orientational order was monitored 
by computing two order parameters. The rank-one order 
parameter Pi is defined as 



Pi = 



Np 
= M 



(37) 



where M|| = M • n is the projection of the total dipole 
moment along the director n (Pi = 1 for a completely 
polarized system). The second-rank order parameter P2 
is the largest eigenvalue of the matrix 



1 



^ 1 



(38) 



where fif is the a component of the vector fJ,^. The cor- 
responding eigenvector n is the director. P2 — ^ when 
all dipoles are oriented parallel to n or — n. 

Table InlHsts our results for these order parameters, as 
well as for the dielectric tensor e = (e||,ej_). In a liquid 
crystal with director n, the latter is determined by the 
following generalization of eq. jS)): 



Npfi 



(39) 



and a similar equation for ej_ in terms of the perpendic- 
ular fluctuations ((M^) - {M i_f ) / N . 



When d is increased above the critical distance dc ~ 
0.63f7, the order parameter P2 jumps from essentially zero 
to about almost one, indicating a first order transition to 
a highly orientationally ordered phase. Fig. II II provides 
snapshots of the simulation cell for d — 0.64(t. It is clear 
from the snapshots that the molecules are associated into 
columns, composed of chains of dipoles oriented head-to- 
tail. These columns are all parallel to the director, and 
are arranged in a hexagonal lattice in the perpendicu- 
lar plane. The simulations for d > O.QAa yielded similar 
liquid crystal phases with columnar order, each with a 
different orientation of the director. The system shows 
strong spatial correFlations in the direction of the direc- 
tor, but it is still fluid in that direction, as indicated by 
the mean-square displacement of the molecules. The lat- 
ter increases indeed linearly with time, with a diffusion 
constant of about _D|| ~ 0.09 • lO^^cm^/s. 

In some runs (not listed in TablellT| , the system formed 
two liquid crystal domains characterized by different ori- 
entations of the director. As the transition to a single 
domain is expected to occur on a time scale much larger 
than our simulation time (8 ns), since it requires the col- 
lective motion of many molecules, we included in Table Hll 
only results from runs where a single domain formed 
spontaneously in less than 10 ns. Most runs yielded fully 
polarized liquid crystals (Pi close to 1). It is likely that 
the lower value of Pi measured in the case d = 0.65(7 
is due to insufficient sampling of phase space: the prob- 
ability of a column inverting its orientation during our 
simulation time is indeed very small. 

As the column configuration is energetically more 
favourable for extended than for point dipoles (see 
Fig- it is not surprising that liquid crystal colum- 
nar phases form at a much lower dipolar coupling con- 
stant than previously reported for point dipolar fluid 
models; here p*"^ /T* = p? /a'^kT ~ 3, while columnar 
phases were observed in the dipolar soft-sphere fluid at 
[i*"^ jT* ~ 9, and in the dipolar hard sphere fluid at 
ji*"^ /T* — 6.25 ^t the same state point as un- 

der consideration here, Kusalik found however that even 
a strong electric field did not induce a nematic phase 
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^ — 









































FIG. 11: Snapshots of the simulation cell of a dipolar soft- 
spheres fluid in the columnar phase. The hexagonal lattice in 
the plane orthogonal to the director is apparent in the first 
snapshot. The dipoles are represented by a line joining the 
minus charge (shown as a small bead) to the plus charge. 



in the dipolar soft-sphere model with point dipoles |35| . 
The hexagonal lattice arrangement found here is to be 
contrasted with the square lattice reported in Ref. Q 
(the lattice type was not specified in Ref. 0). 



VI. CONCLUSION 

We have extended the considerable body of earlier 
work on dipolar fluids by replacing the usual point-dipole 
on spherical molecules by physically more relevant ex- 
tended dipoles obtained by placing two opposite charges 
symmetrically with respect to the centre. The structural, 
dielectric and dynamical behaviour was monitored as the 
spacing d of the charges was increased, keeping the total 
dipole moment fi = qd fixed. Periodic boundary condi- 
tions were used with proper Ewald summations of the 
Coulombic interactions within an infinitely large sphere 
bounded by a dielectric medium of permittivity e'. The 
key findings are the following: 

a) Runs of several nanoseconds, longer than in most 
previous studies, were required to obtain estimates of the 
dielectric constant e within about 2% with the standard 
fluctuation formula l@J. A careful error analysis shows 
that "metallic" boundary conditions (where e' = oo at 
infinity) yield optimal estimates of e. 

b) The values of e deduced from the h^^^ and /i^^" 
correlation functions agree with the fluctuation results 
within statistical errors, provided a sufficiently large sim- 
ulation cell is used to obtain proper estimates of the 
asymptotic behaviour of these correlation functions. The 
strong influence of the boundary condition e' on the pro- 
jection /i^^°(r) arises from the interaction of the polar 



molecules with the reaction fleld due to the dielectric dis- 
continuity between the fluid and the external medium e'. 
When e' 7^ e, /i^^"(r) does not decay to zero at large 
distances, but rather to a finite constant of order 1/V. 
We derived the value of this constant [eq. H22() ] by using 
simple physical arguments based on macroscopic electro- 
statics and linear response theory. 

c) e has a non-trivial dependence on the ratio d* = d/a. 
Up to d* ~ 0.25, e agrees with the point-dipole result 
within statistical uncertainties, thus illustrating the prac- 
tical usefulness of the point dipole limit. For d* > 0.3, 
e drops to a minimum value roughly 6% below the point- 
dipole result when d* ~ 0.55. When d* is further in- 
creased, e increases sharply and reaches a maximum at 
d* ~ 0.6. 

d) For still larger extensions d* , the system is seen 
to undergo a transition, at constant pressure, to an ori- 
entationally ordered stated similar to a columnar phase 
with a hexagonal ordering in the plane orthogonal to the 
director. This phase is characterized by large values of 
the usual orientational order parameters Pi and P2- At 
the same time the dielectric tensor becomes anisotropic, 
and the mean dielectric constant is very low (e ~ 1-4), 
signalling the strong suppression of dipole moment fiuc- 
tuations. The transition to the columnar phase occurs at 
a value of the dipole moment well below that req uired to 
observe the transition with point dipoles [1(tL l34| . 

e) The dynamics, characterized by the relaxation times 
tm and of the total and individual dipole moments, 
as well as by the self diffusion constant D, slows down 
very significantly as d* increases, due to the enhanced 
tendency of the system to form parallel strings, which 
eventually lead to the columnar phase. In the latter, the 
diffusion coefficient parallel to the director is about 
two order of magnitude smaller than D in the isotropic 
phase, but still substantial, showing that the system be- 
haves like a one-dimensional fiuid. 

The present results are for a single high density p* = 
0.8, and a single pressure in the anisotropic phase [p* ~ 
0.22, corresponding to (p*) ~ 1]. Clearly more work is 
needed to be able to map out a complete phase diagram, 
in view of the additional variable d* . The present work il- 
lustrates the strengths and deficiencies of the point-dipole 
model. Many simple molecular systems fall in a range 
d* ~ 0.5, where the deviations from point dipole be- 
haviour begin to be substantial. 
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